Linear operators

Hamiltonian generation

The hamiltonian function can be used to generate a Chern insulator hamiltonian using the following formula:

\[\hat{H} = \sum_i m_i c^\dagger_i \sigma_z c_i + \sum_{x-links} c^\dagger_i \frac{\sigma_z - i \sigma_x}{2} c_j + \sum_{y-links} c^\dagger_i \frac{\sigma_z - i \sigma_y}{2} c_j + h. c.\]

The only essential parameter is m_lattice - it is a CoordinateRepr that defines the $m$ parameter for each site. There is an alternative way to set the m_lattice - with two arguments, a Matrix and a Symbol representation specifier.

Here is an example of both usages:

m_lattice = ones(15, 15)
m_lattice[6:10, 7:11] .= -1
m_repr = CoordinateRepr(m_lattice, :n)
H1 = hamiltonian(m_repr)
H2 = hamiltonian(m_lattice, :n)

H1 == H2 # These are equal ways to specify the m_lattice parameter

true

There is an ability to isolate some sites from other ones. To do that, a domain_mapping matrix is required. It is a CoordinateRepr{Symbol} object, which maps lattice sites to symbols, where different symbols mean different domains, and hoppings between sites in different domains are set to zero. You can define the mapping in the hamiltonian function call, or apply it to an existing hamiltonian matrix using the domains! function:

m_lattice = ones(15, 15)

mapping = CoordinateRepr(fill(:domain1, 15, 15))
mapping[6:10, 6:10] .= :domain2

H1 = hamiltonian(m_lattice, :n, domain_mapping = mapping)
H2 = hamiltonian(m_lattice, :n)
domains!(H2, mapping)

H1 == H2 # These are equal ways to specify the domain mapping

true

You also can apply some magnetic field to the Chern insulator. This can be done by passing a function that takes the position vector and returns the gauge vector potential. You can write your own function use one of the magnetic field macros defined:

• @landau for Landau gauge field
• @symm for symmetric gauge field
• @flux for flux quantum

You can define the field in the hamiltonian function or apply it to an existing hamiltonian matrix using the field! function:

m_lattice = ones(15, 15)
B = 0.1

H0 = hamiltonian(m_lattice, :n, field = (r -> [0, r[1] * B, 0]))
H1 = hamiltonian(m_lattice, :n, field = @landau(B))
H2 = hamiltonian(m_lattice, :n)
field!(H2, @landau(B))

H0 == H1 == H2 # These are equal ways to specify the field

true

Density and coordinate operators

To calculate a density matrix for zero temperature, you can use the filled_projector function. It accepts the hamiltonian matrix and returns a density matrix. You can also specify the Fermi level and the temperature (which are zero by default) - consult the corresponding page at the Library.

You can also generate the coordinate operators using the coord_operators function. It takes an optional parameter symm which defines if the operators are defined symmetrically (i. e. the point with (0, 0) coordinates is in the center of the lattice) or not (i. e. the point with (0, 0) coordinates is in the bottom-left angle of the lattice).

Here is an example of usage of both functions to calculate the local Chern marker:

using Plots, LinearAlgebra

m_lattice = ones(15, 15)
m_lattice[6:10, 7:11] .= -1
H = hamiltonian(m_lattice, :n)

P = filled_projector(H)
X, Y = coord_operators()
c = 4pi * im * P * X * (I - P) * Y * P

heatmap(heatmap_data(c), color = :viridis)

Electric current

You can calculate electric currents using the currents function. It accepts the hamiltonian and the density matrix to produce a matrix with currents. Here is an example of electric currents in a magnetic field:

using Plots

B = 0.1
H = hamiltonian(ones(15, 15), :n, field = @symm(B))
P = filled_projector(H)
ps, qs = quiver_data(currents(H, P) * 10 / B)
quiver(ps, quiver = qs)

Bianca-Resta currents

There are several ways to calculate currents of the local Chern marker. For example, if we define it using the Bianca-Resta formula $c(r) = 4\pi i \langle r | PXQYP | r \rangle$, we can define the current $J(r, r')$ as some formula that complies the following rules:

\begin{cases} dt c(r) = \sum{r'} J(r, r') \ J(r, r') = - J(r', r) \end{cases}, \hspace{0.5cm} d_t c(r) = 4\pi i \langle r | i[H, PXQYP] | r \rangle = 4\pi i \langle r | i([H, P]XQYP + PX[H, Q]YP + PXQY[H, P]) | r \rangle $

Let us call them Bianca-Resta currents.

There are four formulas like that in this package. Each is a macro that takes matrices of the hamiltonian, the density and coordinate operators to generate a function that calculates the current between two sites (given their indices as parameters).

The difference between them is that some of them are translationally invariant - if we shift the coordinate operators (e. g. redefine $\hat{X}$ as $\hat{X} + X_0$), the values of currents will not change. Also some of them are stable - all currents are equal to zero if $P$ is a density matrix of a stationary state of hamiltonian $\hat{H}$ that does not depend on time.

These traits of macros in this package are in the following table:

To find out more about Bianca-Resta currents, check out this section.

Each one takes and returns a function that calculates a current given indices of 2 sites.

To generate a matrix with currents, use the @currents macro.

Here is an example animation of LCM behavior under after hamiltonian quench

using Plots

H1 = hamiltonian(ones(15, 15), :c)
H2 = hamiltonian(ones(15, 15) * 3, :c)

h(t) = H2
P0 = filled_projector(H1)
X, Y = coord_operators()
a = Animation()
@evolution [
    :ham => h => H,
    P0 => h => P
] for t in 0:0.1:3
    c = -4pi * im * P * X * P * Y * P
    cur = @currents @J_treq H P X Y
    plot_auto("LCM currents" => c => cur / 5, hmapclims = (-4, 1))
    frame(a)
end

gif(a, "example_animation.gif", fps = 10)

Streda currents

Another way to calculate the LCM currents is the Streda formula.

This formula defines the LCM as the linear response of the density to the adiabatic field application: $\mathcal{C}_s(r) = \frac{\delta P(r)}{\delta B}$. Such approach allows to evaluate currents of LCM quite easily: $J_{s}(r, r') = \frac{\delta J(r, r')}{\delta B}$.

Note that unlike Bianca-Resta currents, these currents are localized - the current between non-adjacent sites is always zero.

These currents can be evaluated as follows:

τ = 100
time_domain = 0:0.5:2τ

TopologicalMarkers._set_lattice_size!(22, 21)
X, Y = coord_operators()

m0 = fill(-3, 22, 21)
m0[6:12, 7:13] .= -1
m1 = fill(-3, 22, 21)
m1[7:13, 7:13] .= -1
h0 = hamiltonian(m0, :c)
h1 = hamiltonian(m1, :c)
h(t) = h0 + (h1 - h0) * min(1, t/τ)

B = 1e-6
hb0 = hamiltonian(m0, :c, field = @landau(B))
hb1 = hamiltonian(m1, :c, field = @landau(B))
hb(t) = hb0 + (hb1 - hb0) * min(1, t/τ)

P0 = h(0) |> filled_projector
P0b = hb(0) |> filled_projector

a = Animation()
@evolution [
    :ham => h => H,
    :ham => hb => Hb,
    P0 => h => P,
    P0b => hb => Pb,
] for t in time_domain
    cur = currents(H, P)
    curb = currents(Hb, Pb)
    plot_auto("Streda currents" => (Pb - P) / B => (curb - cur) / B,
        hmapclims = (-3, 3), title = "Time: $t")
    frame(a)
end

gif(a, "streda.gif", fps=20)